Australian  Government 
Department  of  Defence 

Defence  Science  and 
Technology  Organisation 


Modelling  of  PBX-115  Using 
Kinetic  CHEETAH  and  the 
DYNA  Codes 

Jing  Ping  Lu  and 
David  L.  Kennedy 

DSTO-TR-1496 


^distribution  Unlimited 


Reproduced  From 
Best  Available  Copy 


20040412  028 


Australian  Government 
Department  of  Defence 

Defence  Science  and 
Technology  Organisation 

Modelling  of  PBXW-115  Using  Kinetic  CHEETAH 
and  the  DYNA  Codes 

Jing  Ping  Lu*  and  David  L.  Kennedy** 

*  Weapons  Systems  Division 

Systems  Sciences  Laboratory 

**  Orica  Explosives 
DSTO-TR-1496 

ABSTRACT 

PBXW-115,  a  highly  non-ideal  explosive  tailored  for  xmderwater  effects  and  composed  of  43% 
ammonium  perchlorate,  25  %  aluminium,  20%  RDX  and  12%  HTPB  binder,  has  been  studied 
using  Kinetic  CHEETAH  (the  Lawrence  Livermore  National  Laboratory  CHEETAH  2.0  code),  the 
three-term  Ignition  and  Growth"  Model  incorporated  into  the  explicit  finite  element  hydrocode 
LSTC-DYNA,  and  the  CPeX  Reactive  Model  based  on  the  small  divergent  detonation  theory  of 
Wood  and  Kirkwood  in  DYNA2D.  This  report  firstly  focuses  on  a  series  of  simulations 
performed  using  Kinetic  CHEETAH  based  on  the  Wood-Kirkwood  detonation  theory  using  a 
pressure  exponent  of  2  in  the  pressure-dependent  rate  law.  It  was  foimd  that  CHEETAH  could 
predict  the  trend  of  the  detonation  velocity  as  a  function  of  charge  diameter,  but  overestimated 
the  detonation  velocities.  The  agreement  was  improved  with  further  parameter  adjustment,  by 
decreasing  the  pressure  exponent  in  the  rate  law  from  2.0  to  1.0,  and  then  0.5.  The  reaction  zone 
widths  over  a  wide  range  of  charge  radius  were  also  computed.  Attention  was  then  turned  to 
hydrocode  modelling,  with  particular  interest  in  developing  reactive  models  for  PBXN-111  and 
PBXW-115  (Aust).  Both  the  "Ignition  and  Growth"  Model  and  the  CPeX  Model  were  calibrated 
against  the  experimentally  observed  dependence  of  detonation  velocity  on  charge  diameter  in 
unconfined  charges  of  both  PBXW-115  (Aust)  and  US  PBXN-111.  These  two  reactive  models  were 
validated  by  comparing  their  predictions  against  experimental  data  for  detonation  of  charges 
confined  in  2.5  mm  and  3  mm  thick  brass.  The  Ignition  and  Growth  Model  was  then  applied 
successfully  to  the  simulation  of  aquarium  tests.  Finally,  to  test  whether  the  two  models  using 
parameters  derived  from  small-scale  tests  can  be  applied  to  large-scale  devices,  these  two  models 
were  applied  to  the  simulation  of  mid-scale  underwater  tests  of  PBXW-115  {Aust),  and 
comparisons  with  available  data  are  made.  The  LS-DYNA  simulations  of  shock  front  curvature 
for  unconfined  charges  are  also  presented  and  the  relationships  between  the  radii  of  curvature  for 
the  detonation  fronts  and  self-propagating  detonation  velocities  are  discussed. 
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Modelling  of  PBXW-115  Using  Kinetic 
CHEETAH  and  the  DYNA  Codes 


Executive  Summary 

In  the  energetic  materials  field,  purely  relying  on  experimental  tests  and  trials  to 
evaluate  the  performance  and  to  vmderstand  the  behaviour  of  new  explosive 
compositions  is  no  longer  an  efficient  approach,  because  they  are  usually  very  costly, 
and  often  unpractical  or  even  impossible  for  complex  configurations.  Therefore,  there 
is  a  continuing  need  for  a  reliable  numerical  modelling  technique  as  an  alternative  tool. 
PBXW-115,  a  extremely  non-ideal  explosive  tailored  for  underwater  effects  and 
composed  of  43%  ammonium  perchlorate,  25  %  aluminium,  20%  RDX  and  12%  HTPB 
binder,  is  a  challenging  explosive  to  model. 

This  report  firstly  focuses  on  a  series  of  simulations  performed  using  Kinetic 
CHEETAH  based  on  the  Wood-Kirkwood  detonation  theory  using  a  pressure  exponent 
of  2  in  the  pressure-dependent  rate  law.  It  was  found  that  CHEETAH  could  predict  the 
trend  of  the  detonation  velocity  as  a  function  of  charge  diameter,  but  overestimated  the 
detonation  velocities.  The  agreement  was  improved  with  further  parameter 
adjustment,  by  decreasing  the  pressure  exponent  in  the  rate  law  from  2.0  to  1.0,  and 
then  0.5.  The  reaction  zone  widths  over  a  wide  range  of  charge  radius  were  also 
computed.  Preliminary  results  using  Kinetic  CHEETAH  indicate  that  the  detonation 
velocity  and  the  critical  diameter  are  sensitive  to  the  assumed  AP  decomposition  rate. 
This  might  suggest  that  the  observed  differences  between  the  US  and  the  Australian 
composition  are  due  to  the  particle  size  of  the  AP  rather  than  the  RDX.  Further 
investigation  of  the  effects  of  the  assumed  AP  decomposition  rate  on  the  detonation 
velocity  as  a  function  of  charge  diameter  will  be  carried  out  in  the  future. 

Attention  was  then  turned  to  hydrocode  modelling,  with  particular  interest  in 
developing  reactive  models  for  PBXN-111  and  PBXW-115  (Aust).  Both  the  "Ignition 
and  Growth"  Model  and  the  CPeX  Model,  although  determined  from  the  detonation 
velocity  versus  diameter  effect  experimental  data  for  tmconfined  charges,  could 
faithfully  model  the  reaction  rate  characteristics  for  experiments  carried  out  in  other 
configurations.  This  was  validated  by  comparing  the  predictions  using  both  models 
against  experimental  data  for  detonation  of  charges  confined  in  2.5  mm  and  3  mm 
thick  brass.  The  Ignition  and  Growth  Model  was  then  appHed  successfully  to  the 
simulation  of  aquarium  tests.  Finally,  to  test  whether  the  two  models  using  parameters 
derived  from  small-scale  tests  can  be  applied  to  large-scale  devices,  these  two  models 
were  applied  to  the  simulation  of  mid-scale  underwater  tests  PBXW-115  (Aust),  and 
comparisons  with  available  data  were  made.  The  LS-DYNA  simulations  of  shock  front 
curvature  for  unconfined  charges  are  presented  and  the  relationships  between  the  radii 
of  curvature  for  the  detonation  fronts  and  self-propagating  detonation  velocities  of 
PBXN-111  are  discussed.  Given  that  there  is  a  lack  of  experimental  wave  curvature 
data  for  PBXW-115(AMSf),  the  simulation  results  presented  herein  provide  some  insight 
into  the  detonation  process  for  this  version  of  the  explosive  composition. 
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1.  Introduction 

Since  the  development  of  PBXW-115‘  (43/25/20/12  AP/Al/RDX/HTPB),  now  known 
as  PBXN-111  in  the  US  after  being  fully  qualified  as  an  underwater  explosive  and 
accepted  for  a  possible  Insensitive  Munition  (IM)  warhead  fill  by  the  US  Navy 
[Anderson  and  Leahy  1985],  research  and  development  have  been  extensively 
conducted  by  many  researchers  for  more  than  a  decade. 

•  Forbes  and  co-workers  at  NSWC  [Forbes  et  al.  1989;  Forbes  et  al.  1992]  have 
applied  different  experimental  techniques  to  obtain  detonation  velocity  and 
critical  diameter  for  both  unconfined  and  confined  charges,  corner  turning 
ability  and  shock  front  curvature.  This  has  provided  a  substantial  data  base  on 
detonation  properties  of  this  non-ideal  explosive,  which  can  be  used  to  calibrate 
hydrodynamic  flow  detonation  models  taking  into  account  ignition  and  growth 
of  reaction  and  multiple  energy  release  rates. 

•  The  CPeX  model  (Commercial  Performance  of  Explosives)  developed  by  Kirby 
and  Leiper  [1985]  was  applied  by  Jones  and  Kennedy  [1991]  to  predict  the 
properties  of  unconfined  PBXW-115  and  foimd  excellent  fit  to  the  experimental 
data.  The  authors  also  made  recommendations  for  future  work  including 
confinement  and  aquarium  test  modelling.  The  same  CPeX  model  was 
implemented  in  the  explicit  finite  element  hydrocode  DYNA2D  to  simulate 
time-dependent  reactive  flow  in  non-ideal  explosives  such  as  PBXN-111  for  a 
variety  of  initiation  and  detonation  tests  [Kennedy  and  Jones  1993],  with  the 
results  generally  in  excellent  agreement  wi^  the  experimental  data. 

•  Miller  [1996a]  developed  a  simplified  ignition  and  growth  (SIG)  model 
consisting  of  only  two  adjustable  parameters,  the  ignition  (I)  and  growth  (G) 
rate  constants  determined  by  iterating  these  variables  in  DYNA2D  hydrocode 
simulations  of  the  failure  diameter  and  the  gap  test  sensitivity  until  the 
experimental  values  were  reproduced.  According  to  the  author,  there  is  a  need 
to  investigate  how  the  pressure  effects  on  reaction  front  radius  of  curvature  and 
/or  the  reaction  zone  thickness  can  be  used  in  further  characterizing  the 
chemical  reaction  rate  kinetics  of  detonation.  He  and  his  co-worker  [Miller 
1996b;  Miller  and  Guirguis  1996]  also  reported  a  reactive  flow  model  developed 
for  detonation  and  combustion  processes  for  non-ideal  metallised  explosives 
that  they  implemented  in  DYNA2D.  The  model  was  based  on  a  coupling  of  the 
global  kinetics  for  the  fast  reactions  that  determine  the  detonation  state  and  the 
slow  metal  combustion  kinetics  that  determine  the  late  time  energy  release  after 
the  CJ  state. 


*  The  terms  PBXW-1 15  and  PBXN-1 1 1  are  used  interchangeably.  PBXW-1 15  was  the  name  used  before 
US  Na\7  qualification.  PBXN-1 1 1  is  the  name  used  post-qualification.  PBXW-1  \5{Aust)  is  the  term  used 
to  identify  the  composition  developed  from  indigenous  Australian  ingredients. 
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•  For  a  number  of  explosives  including  PBXN-111,  Souers  [1997]  developed  a 
model  to  relate  the  size  effects  of  a  cylinder  with  its  average  sonic  reaction  zone 
length.  This  model  assumed  an  external  layer  along  the  edge  of  the  cylinder  out 
of  which  all  energy  was  considered  lost.  Two  major  size  effects  were  identified, 
namely  the  "forward"  effect  (the  reduction  of  the  detonation  velocity  with  the 
decrease  of  the  radius  tmtil  detonation  finally  stops)  and  the  "transverse"  effect 
(the  increase  of  curvature  of  the  detonation  front  with  decreasing  radius).  The 
rate  constants  for  time-dependent  modelling  were  then  estimated  from  reaction 
zone  lengths  based  on  either  the  detonation  front  curvature  or  charge  diameter 
effects  [Souers  and  Garza  1998]. 

•  Souers  [1998]  also  published  a  collection  of  detonation  front  curvature  and  size 
effect  data,  plus  the  calculated  parameters  that  constitute  the  input  to  Kinetic 
CHEETAH.  Howard  et  al.  [1998]  applied  the  renormalised  new  product  library 
NEWCl  including  a  three-phase  carbon  equation  of  state  (EOS)  with  the  more 
reliable  Murnaghan  EOS  for  solids  in  CHEETAH  to  PBXN-111  prediction  and 
achieved  good  agreement  with  experimental  results. 

•  More  recently,  Souers  et  al  [2001]  have  fully  described  and  analysed  the 
cylinder  test  for  deriving  detonation  energies,  including  the  relation  between 
streak  camera  and  Fabry-Perot  interferometer  data.  They  introduced 
CHEETAH  3.0  with  its  new  all-Hugoniot  calibration  producing  the  most 
accurate  detonation  energies  to  date.  Taking  the  composite  explosive  PBXN-111 
as  a  special  problem,  they  found  that  the  Reactive  Flow  code  JWL++  with  one 
fully  reacted  JWL  EOS  failed  to  fit  the  data.  By  constructing  a  2-JWL  form  with 
the  first,  fast  rate  describing  tibe  detonation  velocity  and  the  second,  slow  rate 
the  cylinder  energies,  they  have  successfully  applied  the  JWL++  to  the  non¬ 
ideal  explosive  PBXN-111  that  reacts  as  little  as  25%  in  the  cylinder  test. 

An  "Australianised"  version  of  the  US  Navy  IM  underwater  explosive  PBXN-111, 
based  on  indigenous  ingredients,  was  developed  and  designated  PBXW-115  (Aust) 
[Bocksteiner  and  Billon  1991;  Bocksteiner  et  al.  1994].  Research  by  Bocksteiner  and 
Whelan  first  addressed  binder  and  formulation  studies,  compatibility  and 
sensitiveness,  then  detonation  properties  including  velocity  of  detonation,  critical 
diameter,  shock  sensitivity  and  finally  underwater  explosive  performance  [Bocksteiner 
and  Billon  1991;  Bocksteiner  and  Whelan  1994;  Bocksteiner  1996].  They  found  that 
PBXW-115  (Aust)  had  significantly  different  explosive  properties,  particularly  critical 
diameter  and  shock  sensitivity,  from  its  US  coimterpart,  believed  due  to  the  particle 
size  distribution,  grade  and  type  of  RDX  used  in  the  two  formulations.  These 
investigations  resulted  in  PBXW-115  {Aust)  being  qualified  by  Australian  Ordnance 
Council  for  use  as  an  insensitive  main  charge  fill  in  large  underwater  blast  weapons, 
and  the  proposal  that  the  material  could  be  qualified  as  an  Extremely  Insensitive 
Detonating  Substance  [Whelan  and  Bocksteiner  1998]. 

Turning  attention  to  the  modelling  of  the  Australian  version  of  PBXW-115,  i.e.  PBXW- 
115  (Aust),  the  limited  work  reported  so  far  [Jones  and  Kennedy  1991]  is  confined  only 
to  modelling  detonation  velocity  of  unconfined  PBXW-115  {Aust)  where  the  model  was 
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based  on  the  results  of  experiments  on  relatively  small  charges.  In  order  to  validate 
extrapolation  based  on  the  experimental  results  on  smaller  charges  to  the  larger 
explosive  charges  used  in  imderwater  devices,  the  Weapons  Systems  Division  at 
DSTO-Edinbugh  has  been  undertaking  investigations  into  both  modelling  of  realistic 
charge  geometries  and  small  scale  to  mid-scale  imderwater  tests  to  record  the  "early- 
time"  behaviour  and  performance  of  PBXW-115  {Aust).  The  results  of  the  small-scale 
aquarium  tests  presented  by  Dorsett  and  Katselis  [1999]  and  the  mid-scale  imderwater 
tests  reported  by  Wilkinson  [2001]  are  two  of  the  experimental  investigations  carried 
out,  which  can  be  used  for  calibrating  computer  models.  A  parallel  modelling  program 
was  initially  focussed  on  performance  predictions  with  Lawrence  Livermore  National 
Laboratory  CHEETAH  2.0  code  based  on  both  the  traditional  Chapman-Jouget  (CJ) 
thermod5mamic  detonation  theory  and  the  Wood-Kirkwood  kinetic  detonation  theory 
[Lu  2001]. 

Following  the  earlier  studies  of  non-ideal  explosive  PBXW-115  performance  [Lu  2001; 
Lu  and  Kermedy  2001],  and  in  order  to  further  characterize  the  non-ideal  features  of 
this  explosive,  this  report  firstly  focused  on  a  series  of  simulations  performed  using 
Kinetic  CHEETAH  based  on  the  Wood-Kirkwood  detonation  theory  using  a  pressure 
exponent  of  2  in  the  pressure-dependent  rate  law.  It  was  found  that  CHEETAH  could 
predict  the  trend  of  the  detonation  velocity  as  a  function  of  charge  diameter,  but 
overestimated  the  detonation  velocities.  The  agreement  was  improved  with  further 
parameter  adjustment  by  decreasing  the  pressure  exponent  in  the  rate  law  from  2.0  to  1 
and  then  0.5.  The  reaction  zone  widths  over  a  wide  range  of  charge  radius  were  also 
computed. 

Attention  was  then  turned  to  hydrocode  modelling,  with  particular  interest  in 
developing  reactive  models  for  PBXN-111  and  PBXW-115  {Aust)  [Lu  et  al  2002].  In  this 
report,  both  the  three-term  "Ignition  and  Growth"  Model  incorporated  into  the  explicit 
finite  element  hydrocode  LS-DYNA,  and  the  CPeX  Reactive  Model  based  on  the  small 
divergent  detonation  theory  of  Wood  and  Kirkwood  in  DYNA2D,  have  been  calibrated 
against  the  experimentally  observed  dependence  of  detonation  velocity  on  charge 
diameter  in  imconfined  charges  of  both  PBXW-115  (Aust)  and  US  PBXN-111.  These  two 
reactive  models  were  validated  by  comparing  their  predictions  with  experimental  data 
for  detonation  of  charges  confined  in  2.5  mm  and  3  mm  thick  brass.  The  Ignition  and 
Growth  Model  was  then  applied  successfully  to  simulation  of  aquarium  tests.  Finally, 
to  test  whether  the  two  models  using  parameters  derived  from  small-scale  tests  can  be 
applied  to  large-scale  devices,  these  two  models  were  applied  to  the  simulation  of  mid¬ 
scale  underwater  tests  on  PBXW-115  (Aust),  and  comparisons  with  available  data  were 
made.  The  LS-DYNA  simulations  of  shock  front  curvature  for  unconfined  charges  have 
also  been  presented  and  the  relationship  between  the  radii  of  curvature  for  the 
detonation  fronts  and  self-propagating  detonation  velocities  has  been  discussed. 
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2.  Kinetic  Calculations 


Kinetic  CHEETAH  is  based  on  the  Wood-Kirkwood  (WK)  detonation  theory  [Wood 
and  Kirkwood  1954]  which  is  specially  designed  for  modelling  time-dependent 
phenomenon.  The  new  chemical  kinetics  model  implemented  in  CHEETAH  considers 
detonation  in  composite  explosives  with  large  reaction  zones,  and  the  interplay 
between  the  energy  produced  by  kinetically  controlled  reactions  and  the  energy  lost 
due  to  radial  expansion  of  the  product  gases.  Wood-Kirkwood  theory  thus  allows 
prediction  of  the  dependence  of  detonation  parameters  on  charge  diameters,  and 
estimation  of  the  length  of  the  detonation  zone,  identified  as  the  region  behind  the 
detonation  wave  for  which  the  sum  of  the  mass  velocity  and  the  velocity  of  sormd  is 
equal  to  the  detonation  velocity  [Loboiko  and  Lubyatinsky  2000]. 

As  described  in  the  CHEETAH  2.0  User's  Manual  [Fried  et  al  1998],  WK  theory  starts 
with  the  hydrodynamic  Euler  equations  coupled  to  chemical  kinetics.  The  theory  treats 
the  detonation  along  the  centre  of  the  cylinder.  Radial  expansion  is  treated  as  a  first 
order  perturbation  to  perfect  one  dimensional  planar  detonation.  The  Euler  equations 
are  reduced  to  their  steady  state  form.  The  result  is  a  set  of  ordinary  differential 
equations  that  describe  hydrodynamic  variables  and  chemical  concentrations  along  the 
centre  of  the  cylinder.  The  theory  requires  specification  of  the  rate  of  radial  expansion, 
(Or,  as  a  function  of  radius.  Although  Kinetic  CHEETAH  has  implemented  three  radial 
expansion  models  in  the  code,  in  this  study  the  simple  pressure  model  with  the 
following  time  rate  of  change  of  (Or  is  used: 

d  0)r  2SP  p  2 

—J— - 

dt  Rip  a  (1) 

where  cOr  (f  =  O)  =  (Z>^  -u)lR^  (2) 

Here,  P  is  the  pressure,  u  is  the  particle  velocity  in  the  shock  frame,  po  is  the  initial 
density  of  the  explosive,  Ro  is  the  charge  radius  and  S  is  an  empirical  scaling  factor.  If 
this  model  is  used  with  S  =  0,  is  a  constant  with  the  initial  value  determined  by  the 
radius  of  curvature  Rc,  the  detonation  velocity  Ds,  and  the  particle  velocity  at  the 
detonation  front.  The  radius  of  curvature  is  obtained  from  Souer's  detonation  front 
curvature  and  size  effect  data  [Souers  1998]. 

Kinetic  CHEETAH  assumes  the  concentrations  of  individual  reactants  are  controlled 
by  the  rate  of  the  kinetic  reactions,  while  the  products  are  assumed  to  be  in 
thermochemical  equilibrium.  Kinetic  CHEETAH  supports  multiple  reaction  rate  laws: 

•  Simple  constant  reaction  rate  law 

•  Simple  Arrhenius  kinetics  with  a  temperature-dependent  pre-factor 

•  Pressure-dependent  rate  law 

•  Hot  spot  model 
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A  simple  constant  reaction  rate  has  the  form: 

%=Ra-X)  (3) 

where  R  is  the  rate  constant,  and  A  represents  the  amount  of  unbumed  reactant 
normalised  to  vary  between  zero  (all  tmburned)  and  one  (all  burned). 

The  Arrhenius  decay  law  is  : 

%  =  «r’'exp(-B/r)(l-A))  (4) 

where  T  is  the  temperature,  R  and  B  are  the  rate  constants,  and  A  is  the  temperature 
exponent. 

A  pressure-dependent  rate  law  describing  surface  controlled  reaction  has  the  form: 

=  (5) 

where  P  is  the  pressure,  R  is  the  rate  constant,  A  and  B  are  the  extent  of  reaction 
exponents,  and  D  is  the  pressure  exponent. 

The  hot  spot  model  has  the  form: 

where  E  and  C  are  constants,  and/ represents  the  degree  of  reaction  of  hot  spots 

CHEETAH  has  used  the  following  simple  pressure-dependent  rate  law  to  infer  kinetic 
rates  for  a  variety  of  high  explosives  and  their  composites: 

‘‘y^,=(}-x)RP^  (8) 

We  have  also  used  the  same  rate  law  in  our  study.  By  using  the  NEWCl  product 
library  and  the  updated  rate  constant  values  for  the  pressure  exponent  D  of  2  as 
defined  by  Howard  et  al.  [1998],  rather  than  those  rate  constants  initially  recorded  in 
the  CHEETAH  2.0  User's  Manual,  we  have  computed  the  reaction  zone  widths  over  a 
wide  range  of  charge  radius,  and  these  results  are  presented  in  Table  1.  Souers  [1999] 
proposed  a  dimensionless  relation  N  ^xjR^to  measure  the  non-ideality  of  explosives 

(where  Xe  is  the  reaction  zone).  For  an  ideal  explosive,  N  =  0.  The  larger  N  gets,  the 
more  non-ideal  the  explosive  becomes.  For  N  >  0.3,  the  explosive  is  extremely  non¬ 
ideal.  According  to  this  definition,  PBXW-115  can  be  classified  as  extremely  non-ideal 
for  Ro  <  24.8  mm,  moderately  non-ideal  for  24.8  mm  <  Ro  <  85  mm  and  mildly  non¬ 
ideal  for  Ro  >  100  mm  (see  Table  1). 


(6) 

(7) 
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Table  1.  Reaction  zones  calculated  from  Kinetic  CHEETAH 


Charge  radius 
Ro  (mm) 

Time  to  sonic  point 
(lis) 

Reaction  zone 
;Ce(mm) 

N  =  xJR, 

19.35 

1.315 

6.786 

0.351 

20.55 

1.288 

6.814 

0.332 

24.8 

1.279 

7.072 

0.285 

34.55 

1.291 

7.410 

0.214 

49.95 

1.347 

7.917 

0.159 

65 

1.370 

8.173 

0.126 

75 

1.361 

8.160 

0.109 

85 

1.455 

8.759 

0.103 

100 

1.498 

9.054 

0.091 

Using  the  same  rate  constants  with  the  pressure  exponent  D=2,  we  also  calculated  the 
detonation  velocities  as  a  function  of  reciprocal  of  radius  of  curvature.  Figure  1 
summarises  the  predicted  results  with  Kinetic  CHEETAH,  compared  with  the 
experimental  data  [Forbes  et  al  1989;  1992  and  Souers  1998].  It  is  seen  that  although 
CHEETAH  can  not  reproduce  exactly  the  experimental  results,  it  follows  a  similar 
trend  which  demonstrates  the  diameter  effect.  The  agreement  could  be  improved  with 
further  parameter  adjustment. 


Figure  1  Detonation  velocihj  versus  reciprocal  radius  of  curvature  ofPBXW-115  with  D=2 


By  calibrating  kinetic  parameters  based  on  the  detonation  velocity  versus  diameter 
effect  experimental  data  for  unconfined  charges  [Forbes  et  al.  1989;  Forbes  et  al.  1992], 
we  have  developed  the  rate  constants  by  varying  the  pressure  exponent  D  in  the  rate 
law  from  2.0  to  1  and  0.5.  The  rate  constants  R  for  the  reactants  developed  in  this  study 
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are  listed  in  Table  2.  The  predicted  detonation  velocities  with  the  pressure  exponent  D 
of  1  and  0.5  together  with  the  experimental  data  are  given  in  figures  2  and  3.  It  is  seen 
that  the  agreement  was  much  improved  by  decreasing  the  pressure  exponent  D  in  the 
rate  law  from  2.0  to  0.5.  Northam  and  Jessee  [1969]  conducted  experimental 
investigations  of  the  effect  of  aluminium  size  and  loading  on  die  burning  rate  of  solid 
propellants  (AP/  Al/ curing  agents)  under  various  acceleration  loads.  At  each 
acceleration  level  the  average  burning  rates  were  fitted  by  the  least-squares  method  to 
Vieille's  burning  rate  law  [r=a  (P/500)^  ,  where  r  is  the  average  burning  rate  over 
pressure  action  time  in  mm/s,  n  is  the  pressure  exponent,  P  is  the  average  chamber 
pressure  in  N/ m2  and  a  is  a  constant).  The  pressure  exponent  n  resulted  from  the  data 
fit  ranged  from  0.239  to  0.491.  Their  experimental  data  lend  some  support  to  the 
square  root  pressure  dependence  providing  a  better  fit  to  the  detonation  velocity  data 
than  the  pressure  exponent  of  either  1.0  or  2.0. 


Table  2.  Rate  constant  R  used  in  pressure-dependent  rate  laws 


Reactant 

R  (|u.s-iGPa-2) 

D=2 

D=1 

D=0.5 

Al 

0.0075 

0.084375 

0.3 

AP 

0.0075 

0.084375 

0.3 

HTPB 

0.001 

0.01125 

0.04 

RDX 

0.2 

2.25 

12.5 

Figure  2  Detonation  velocity  versus  reciprocal  radius  of  curvature  ofPBXW-115  with  D=1 
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Figure  3  Detonation  velocih/  versus  reciprocal  radius  of  curvature  ofPBXW-115  ivith  D=0.5 


Preliminary  modelling  results  conducted  by  the  second  author  indicate  that  the 
detonation  velocity  and  the  critical  diameter  are  sensitive  to  the  assumed  AP 
decomposition  rate.  Although  this  might  suggest  that  die  observed  differences  between 
the  US  and  the  Australian  composition  are  due  to  the  particle  size  of  the  AP  rather  than 
the  RDX,  it  is  worth  noting  that  the  recent  study  on  the  effect  of  Australian  RDX  on 
PBXN-109  performances  [Lochert  et  al  2002]  indicated  that  the  above  differences  (large 
critical  diameter  and  improved  shock  sensitivity)  were  also  present  in  PBXN-109 
(which  has  no  AP).  To  clarify  the  role  of  AP,  further  investigation  of  the  effects  of  the 
assumed  AP  decomposition  rate  on  the  detonation  velocity  as  a  fxmction  of  charge 
diameter  will  be  required  in  the  future. 


3.  Hydrocode  Modelling 

3.1  Ignition  and  growth  reactive  model  in  LS-DYNA 

The  Ignition  and  Growth  (I  &  G)  Reactive  Model  in  LS-DYNA[1999]  is  derived  from 
the  original  work  of  Lee  and  Tarver  [1980],  which  permits  the  resolved  reaction  zone 
simulation  of  initiation  (or  failure  to  initiate)  and  detonation  wave  propagation  of  solid 
high  explosives.  According  to  Tarver  et  al.  [1985],  shock  initiation  of  heterogeneous 
solid  explosives  should  be  modelled  as  at  least  a  three-step  process.  The  first  step  is  the 
formation  of  hot  spots  created  by  various  mechanisms  (void  closure,  viscous  heating, 
shear  banding,  etc.)  during  shock  compression  and  the  subsequent  ignition  (or  failure 
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to  ignite  due  to  heat  conduction  losses)  of  these  heated  regions.  The  second  step  in  the 
process  is  assumed  to  be  a  relatively  slow  growth  of  reaction  in  inward  and/or 
outward  "burning"  of  the  isolated  hot  spots.  The  third  step  in  the  shock  initiation 
process  is  a  rapid  completion  of  the  reaction  as  the  reacting  hot  spots  begin  to  coalesce. 
This  model  requires  [Murphy  et  al.  1993]: 

•  An  unreacted  explosive  equation  of  state; 

•  A  reaction  product  equation  of  state; 

•  A  reaction  rate  law  that  governs  the  chemical  conversion  of  explosive 
molecules  to  reaction  product  molecules; 

•  A  set  of  mixture  equations  to  describe  the  states  attained  as  the  reactions 
proceed. 

Both  imreacted  and  product  equations  of  state  are  of  Jones-Wilkins-Lee  (JWL)  forms 

where  P  is  pressure,  V  is  relative  volume,  T  is  temperature,  and  A,  B,  Ri,  R2,  co  (the 
Grxmeisen  coefficient)  and  Cv  (the  average  heat  capacity)  are  constants.  The  chemical 
reaction  rate  equation  in  the  three-term  ignition  and  growth  model  is  of  the  form 
[Tarver  and  Green  1989]: 


dF 

dt 


'  P  , 

— — \-a 


Po 


+  Gi(^-FyF^py+G2{i-FyFfp^  (lO) 


where  P  is  the  fraction  reacted,  t  is  time,  po  is  initial  density,  p  is  current  density,  P  is 
pressure,  and  1,  Gi,  G2,  b,  x,  a,  b,  c,  d,  y,  e,f,  and  z  are  constants.  Three  more  constants  are 
added  to  the  model:  Fmxig,  FmxCr  and  FmnGr  which  limit  the  contributions  of  the  three 
terms  to  respectively  a  maximum  reacted  fraction  Fmixg  for  the  first  term,  a  maximum 
fraction  FmxCr  for  the  second  term  and  a  minimum  fraction  FmnCr  for  the  last  term. 


It  is  known  that  the  mesh  size  significantly  influences  the  prediction  accuracy.  The 
mesh  size  dependency  of  the  model  with  different  mesh  densities  were  examined. 
Considering  both  the  factors  of  computing  time  and  accuracy,  an  element  size  of  1mm 
for  explosives  and  2mm  for  water  was  required  to  model  the  experiment  correctly, 
giving  close  to  8  mesh  points  tiirough  the  reaction  zone. 

Table  3  lists  the  parameters  in  the  Ignition  and  Growth  Reactive  Model  fitted  to 
experimental  data  for  the  detonation  velocity  versus  diameter  for  unconfined  PBXN- 
111  [Forbes  1989]  and  PBXW-115  (Aust)  [Bocksteiner  et  al.  1994].  Figure  4  summarises 
the  LS-DYNA  predictions  and  kinetic  CHEETAH  predictions  of  detonation  velocity  in 
axisymmetric  geometry,  together  with  experimental  data.  Agreement  is  seen  to  be 
excellent  between  the  DYNA  predictions  and  the  experimental  results  for  both  PBXW- 
115  (Aust)  and  PBXN-111  (US).  The  curve  predicted  by  LS-DYNA  and  that  predicted  by 
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CHEETAH  follow  similar  trends.  As  mentioned  earlier,  even  though  CHEETAH 
overestimates  the  detonation  velocities,  it  can  still  predict  the  trend  of  the  detonation 
velocity  versus  diameter  effect. 


Table  3.  Parameters  for  the  Ignition  and  Growth  of  Reaction  Model  [Lu  et  al  2002] 


1  Explosive 

PBXN-111 

PBXW-115(Aust)  1 

1  Unreacted  Equation  of  State  and  Constitutive  Values  | 

po  (g/cm3) 

1.792 

A  (GPa) 

4.066E+03 

B  (GPa) 

-133.9 

Ri 

7.2 

Ri 

3.6 

R3=0}*C4GPa/K) 

2.091E-03 

Yield  Strength  (GPa) 

0.2 

Shear  Modulus  (GPa) 

4.54 

1  Reacted  Product  Equation  of  State  and  CJ  Values  | 

A  (GPa) 

372.9 

B  (GPa) 

5.412 

Ri 

4.453 

Ri 

1.102 

R4=0)-rCv  {GPa/K) 

4.884E-04 

Eo(kJ/cc) 

12.95 

Dcj(mm/ps) 

6.476 

Pci  (GPa) 

20.84 

1  Reaction  Rate  Parameters  for  3  Term  Model  | 

I  (psec-i) 

30 

15 

b 

0.6667 

0.6667 

a 

0 

0 

X 

4.0 

4.0 

Gj  (GPa-ypsec'i) 

0.045 

0.0195 

c 

0.6667 

0.6667 

d 

0.1111 

0.1111 

y 

1.0 

1.0 

G2  (GPa-z|xsec-i) 

1.805E-03 

8.00E-04 

e 

1.0 

1.0 

f 

0.1111 

0.1111 

z 

2.0 

2.0 

P  mixg 

0.015 

0.015 

PmxGr 

0.25 

0.25 

PmnGr 

0 

0 

Note:  Eo  is  the  internal  energy. 
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6500 

6300 

6100 

5900 

5700 

5500 

5300 

5100 

4900 

4700 

4500 

0  0.01  0.02  0.03  0.04  0.05  0.06  0.07 

1/r  (1/mm) 

Figure  4  PBXW-115  Detonation  velocity  for  unconfined  charges 

3.2  CPeX  reaction  model 

The  CPeX  (Coirtmerical  Performance  of  eXplosives)  Reaction  Model  is  based  upon  the 
Wood-Kirkwood  theory  of  slightly  divergent  detonation.  Like  the  Ignition  and  Growth 
model,  this  model  also  involves  a  three-stage  process  to  describe  varying  rates  of 
combustion  among  the  ingredients  of  a  composite  explosive.  It  describes  the  flow  along 
the  central  streamtube  between  the  detonation  front  and  the  CJ  plane  for  unconfined 
cylindrical  geometry.  There  are  four  adjustable  parameters  -  three  characteristic 
reaction  times  (hotspot,  intermediate  and  final  stages  of  the  reaction)  and  the  critical 
pressure  that  mhibits  the  onset  of  the  hotspot  reaction,  which  can  be  calibrated  against 
the  experimentally  observed  dependence  of  detonation  velocity  on  charge  diameter. 
For  both  PBXN-111  and  PBXW-115(Ausf),  these  parameters  are  based  on  the 
assumption  of  the  initial  ignition  and  consumption  of  the  RDX,  the  intermediate 
decomposition  of  the  AP  plus  binder,  and  the  later  reaction  of  the  Al.  The  CPeX  model 
has  been  incorporated  into  both  the  explicit  finite  difference  two-dimensional  in-house 
multi-material  Eulerian  hydrocode  "MULTI"  [Jones  et  al  1998]  and  the  explicit  finite 
element  hydrocode  DYNA2D  [Kennedy  and  Jones  1993]  to  permit  the  simulation  of 
time-dependent  reactive  flow  in  non-ideal  explosives. 
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3.3  Simulations  of  detonation  front  curvature 

The  detonation  velocities  of  highly  non-ideal  explosives  can  show  large  departures 
from  the  ideal  detonation  velocities  predicted  from  thermodynamic  equilibrium  codes 
(Kennedy  1998).  This  is  due  to  detonation  waves  in  non-ideal  composite  explosives 
having  curved  wave  fronts  resulting  from  the  reacting  explosive  having  a  finite  size,  a 
varying  energy  release  rate,  and  hydrodynamic  flow  [Forbes  et  al  1992].  Knowledge  of 
the  variation  of  shock  front  curvature  with  charge  diameter  is  important  in  that  it  is 
basic  data  for  validation  and  use  of  detonation  models  [Leiper  and  Hacked  1997]. 
Experimental  wave  curvature  data  for  PBXW-115(AMSf)  is  not  yet  available;  herein  the 
LS-DYNA  simulations  of  shock  front  curvature  for  both  PBXN-111  and  PBXW- 
115(Awsf)  unconfined  charges  using  Ignition  and  Growth  Reactive  Model  are 
presented,  where  the  simulated  results  for  PBXN-111  are  compared  with  the  available 
experimental  data  [Forbes  et  al  1992;  Lemar  and  Forbes  1994].  Analogy  between  the 
two  versions  of  this  explosive  is  expected. 

Although  the  detonation  front  locus  can  be  fitted  to  many  forms,  such  as  circles,  second 
order  polynomials,  ellipses  and  the  natural  logarithm  of  a  Bessel  function  [Kennedy 
1998],  herein  two  forms  are  chosen  to  fit  the  simulated  locus.  The  first  is  an  ellipse  in 
the  form  [Kermedy  1998] 


1-Jl- 


(11) 


where  R  is  the  radial  coordinate  along  the  charge  radius  (see  Figure  5),  L  the  detonation 
lag  along  the  axis  of  the  cylindrical  charge  (see  Figure  5),  and  a  and  b  are  fitting 
constants  for  a  given  radius.  The  second  is  the  equation  described  by  Souers  and  Raul 
[1997] 


L(R)  =  aR2  +  bR^ 


(12) 


The  fitted  results  of  the  simulations  for  PBXN-111  together  with  the  experimental 
results  are  shown  in  Figure  6.  The  scatter  indicates  the  difficulty  in  measuring  the  wave 
front  curvature  experimentally.  It  is  seen  that  the  curves  simulated  using  both 
equations  (11)  and  (12)  reproduce  the  experimental  data  very  well.  As  noted  by  Miller 
and  Sutherland  [1995],  the  calculated  effects  at  the  edges  of  the  charges  on  the 
detonation  front  depend  greatly  on  the  grid  size  used  in  the  hydrocode  computations. 

Curvature  of  wave  front  at  any  radial  coordinate  R  is  defined  by 
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The  radius  of  curvature  is  given  by 

Rc(R)=1/k{R)  (u) 

The  fit  to  equations  (11)  and  (12)  allows  the  radius  of  curvature  at  any  point  on  the 
wave  front  to  be  computed  by  Eqs.  (13)  and  (14).  The  resulting  functions  for  the  three 
PBXN-111  charges  are  plotted  in  Figure  7.  Similar  to  the  Bessel  function  fit  [Forbes  and 
Lemar  1998],  the  magnitude  of  the  radius  of  curvature  reaches  the  maximum  at  the 


Figure  5  A  sketch  diagram  of  a  typical  rate  stick  which  shows  variables  L,  R,  Band  a 


Radial  Coordinate,  R  (mm) 

Figure  6  Detonation  curvature  of  PBXN-111 
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Figure  7  Radius  of  curvature  as  a  function  of  radial  position  for  PBXN-111 

The  angle  between  the  normal  to  the  locus  and  the  axis  (see  Figure  5),  6(R),  is  given  by 

6(R)=  arctan(L'(R))  (15) 

The  angle  between  the  detonation  wave  and  the  edge  of  the  charge,  a,  can  be 
determined  by  two  methods.  In  the  first  method  [Forbes  et  al  1992],  the  simulated  data 
over  the  last  3mm  at  each  edge  of  the  charge  are  least  squares  fitted  to  a  straight  line  to 
obtain  the  angles.  In  the  second  method,  the  derivative  of  either  equation  (11)  or 
equation  (12)  is  used  to  obtain  the  slope  of  the  wave  front  as  given  in  equation  (15), 
from  this,  and  the  angle  a  (see  Figure  5)  is  calculated  by 

a=9(P-d(Ro)  (16) 

where  Ro  is  the  radius  of  the  explosive  charge. 

The  angles  at  the  edges  of  the  charges  with  different  radii,  calculated  by  both  methods, 
and  the  radii  of  curvature  at  the  centre  of  the  charges,  are  summarized  in  Table  4.  For 
comparison,  the  experimental  data  [Forbes  et  al  1992;  Lemar  and  Forbes  1994]  are  also 
included.  It  is  seen  that  the  angles  obtained  from  a  linear  fit  achieve  the  best  correlation 
with  the  experimental  results.  In  general,  the  radii  of  curvature  at  the  charge  centre 
obtained  from  an  ellipse  fit  as  given  in  equation  (11)  results  in  a  better  correlation  with 
the  experimental  data  than  equation  (12),  whereas  angles  obtained  from  equation  (12) 
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fit  give  better  correlation  than  those  obtained  from  an  ellipse  fit.  A  detailed  study  on 
multi-valued  normal  shock  velocity  versus  curvature  relationships  for  highly  non-ideal 
explosives  conducted  by  Kennedy  [1998]  showed  that  it  was  not  straightforward  to 
derive  the  curvature  close  to  the  edges  of  the  charges  and  the  deduced  curvature  was 
strongly  dependent  on  the  details  of  the  fitting  procedure.  Although  we  have 
attempted  to  fit  the  simulated  results  to  a  modified  elliptical  form  by  a  higher  order 
correction  term  as  given  in  [Kennedy  1998],  the  correlation  has  not  been  improved. 


Table  4.  Angles  of  the  detonation  wave  at  the  edge  of  the  charges  and  the  radii  of  curvature  at 
the  wave  centre 


Material 

Diameter 

(mm) 

a+ 

(Deg) 

a* 

(Deg) 

a*' 

(Deg) 

a* 

(Deg) 

Rc* 

(mm) 

Rc** 

(mm) 

Rc## 

(mm) 

PBXN- 

111 

41.1 

58.9 

50.4 

53.6 

59.1 

58.2 

55.9 

60.65 

48.02 

60.5 

52.2 

57.9 

56.9 

78.4 

76 

66.05-73.4 

68.25 

60.9 

50.8 

55.5 

60.7 

126.4 

123 

142.6-149.5 

-I-  Obtained  from  a  linear  fit  of  the  simulated  data  at  the  edge  over  the  last  3mm. 

*  Obtained  from  equation  (11)  fit. 

**  Obtained  from  equation  (12)  fit. 

#  Data  from  [Lemar  and  Forbes  1994]. 

##  Data  from  [Forbes  et  al  1992]  (Rc  =3.18d  —  76.30,  where  d  is  the  charge  diameter). 


Souers  and  Garze  [1998]  presented  a  simple  wave  curvature  theory  to  relate  curvature 
with  the  size  effect  where  they  proposed  the  following  formula  for  calculating  the 
average  sonic  reaction  zone  length,  : 


_Ro(l-Us^  Poo) 
sm(^)cos(^) 

where  is  charge  radius,  IJ ^  is  the  detonation  velocity  and  is  the  infinite 
diameter  detonation  velocity. 

Using  equation  (17),  the  zone  lengths  based  on  the  angles  determined  with  equation 
(12)  are  calculated  for  PBXN-111  and  the  results  are  shown  in  Figure  8.  The  zone 
lengths  are  normalized  with  respect  to  charge  radius /?q.  For  comparison,  the  zone 
lengths  calculated  using  the  angles  obtained  from  experimental  data  [Lemar  and 
Forbes  1994]  for  PBXN-111  and  those  calculated  from  Kinetic  CHEETAH  are  also 
included.  The  agreement  is  reasonable  between  the  DYNA  predictions  and  the 
experimental  results.  Even  though  CHEETAH  underestimates  the  zone  lengths,  the 
curve  predicted  by  LS-DYNA  and  that  predicted  by  CHEETAH  follow  similar  trends. 
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Following  the  above  fitting  procedure,  the  simulated  curves  of  detonation  curvature 
and  radii  of  curvature  as  a  function  of  radial  distances  for  PBXW-115  (Aust)  are  plotted 
in  Figure  9  and  Figure  10,  respectively.  The  angles  at  the  charge  edges  and  the  radii  of 
curvature  at  the  charge  centres  are  given  in  Table  5.  Figure  11  shows  the  detonation 
velocity  as  a  function  of  curvature  at  wave  front  centres.  It  is  seen  that  the  detonation 
velocity  is  sensitive  to  the  curvature  of  the  detonation  front.  The  smaller  the  charge 
diameters,  the  greater  expansion  and  wave  curvature,  hence  result  in  the  lower  the 
detonation  velocity.  Figure  12  presents  the  normalized  average  sonic  reaction  zone 
lengths  calculated  with  equation  (17). 

Given  that  there  is  a  lack  of  experimental  wave  curvature  data  for  PBXW-115(AMsf),  the 
simulation  results  presented  herein  could  provide  some  insight  into  the  detonation 
process  for  this  non-ideal  explosive. 
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Figure  9  Detonation  curvature  ofPBXW-115  (Aust) 


Figure  10  Radius  of  curvature  as  a  function  of  radial  position  for  PBXW-115  (Aust) 
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Table  5.  Angles  of  the  detonation  wave  at  the  charge  edges  and  the  radii  of  curvature  at  the 
wave  centre  for  PBXW-115(Aust) 


Material 

Diameter 

(mm) 

«• 

(Deg) 

a” 

(Deg) 

Rc* 

(mm) 

Rc" 

(mm) 

PBXW- 

115(Awsf) 

80 

51.1 

51.5 

88.4 

84.8 

100 

54.3 

58.2 

173.9 

167 

140 

59.3 

61.4 

245 

236 

170 

59.3 

61.4 

306 

293 

200 

62.3 

62.6 

327.9 

313.4 

*  Obtained  from  equation  (11)  fit.  **  Obtained  from  equation  (12)  fit. 


Figure  11  Detonation  velocity  of  PBXW-115  (Aust)  as  a  function  of  curvature  at  wavefront 
centres 
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Figure  12  Average  sonic  reaction  zone  lengths  ofPBXW-115  (Aust) 

3.4  Modelling  of  confined  charges 

Simulating  confined  charges  is  difficult.  The  detonation  velocity  is  sensitive  to  the 
curvature  of  the  detonation  front.  The  conditions  at  the  explosive/confinement 
interface  have  an  influence  on  the  curvature.  The  element  aspect  ratio  in  the 
confinement  should  be  close  to  1  as  greater  aspect  ratios  lead  to  the  shock  in  the 
confinement  travelling  faster  than  it  should,  and  giving  too  flat  a  detonation  front  in 
the  explosive.  In  order  to  consider  two  surfaces  in  contact  it  is  necessary  to  designate 
one  as  a  slave  surface  and  the  other  as  a  master  surface.  Parts  in  DYNA  keyword  file 
are  composed  of  material  and  section,  equation  of  state  and  hourglass  data.  Parts  in  the 
slave  part  set  are  checked  for  contact  with  parts  in  the  master  part  set.  Self  contact  is 
checked  for  any  part  in  both  sets.  The  Ignition  and  Growth  Reactive  Model  in  LS- 
DYNA  and  the  CPeX  Reactive  Model  in  DYNA2D  as  described  previously  were 
validated  by  comparing  their  predictions  against  experimental  data  [Forbes  1989; 
Bocksteiner  et  al.  1994]  for  detonation  of  charges  confined  in  2.5  mm  and  3  mm  thick 
brass.  Figure  13  demonstrates  the  I  &  G  model  and  the  CPeX  model  predictions  of 
detonation  velocity  in  axisymmetric  geometry,  compared  favourably  with  the 
experimental  data  (for  confined  detonations  in  both  2.5mm  and  3  mm  thick  brass). 
Excellent  agreement  is  also  seen  between  the  Ignition  and  Growth  Model  and  the  CPeX 
Model. 
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Figure  13  PBXW-115  Detonation  velocity  for  confined  charges 

3.5  Simulation  of  aquarium  tests 

The  Ignition  and  Growth  Reactive  Model  described  above  has  been  used  to  simulate 
the  small-scale  aquarium  tests  reported  by  Dorsett  and  Katselis  [1999].  The  simulation 
used  some  63000  elements  in  axisymmetric  cylindrical  geometry  with  a  cell  size  of 
1mm.  The  PE4  booster  was  detonated  using  the  programmed  bum  option  in  LS-DYNA 
[1999]  with  the  JWL  parameters  [Ly  2000]  given  in  Table  6.  The  PBXW-115  (Aust) 
charge  with  a  diameter  of  100mm  and  a  length  of  200mm  was  burnt  by  the  above 
Ignition  and  Growth  Reactive  Model.  The  surrounding  water  was  modelled  with 
Gruneisen  equation  of  state  with  the  parameters  as  given  in  Table  7  [Meyers  1994]. 


Table  6.  JWL  parameters  for  PE4  booster  [Ly  2000] 


po(g/cm3) 

1.59 

A  (GPa) 

774.054 

B  (GPa) 

8.677 

Ri 

4.837 

Ri 

1.074 

Q) 

0.284 

Eo(KJ/cc) 

9.381 

Dq  (mm/ps) 

7.9 

Pcj(GPa) 

24.0 
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Table  7.  Gruneisen  parameters  for  water  [Meyers  1994] 


P)(g/cm3) 

1.00 

To 

0.1 

C{mm/ps) 

1.65 

S 

1.92 

Note:  /o  is  the  Gruneisen  gamma,  C  is  the  intercept  of  the  shock  velocity 
vs  particle  velocity  curve  and  S  is  the  coefficient  of  the  slope  of  this  curve . 

By  tracking  the  initial  high-pressure  shock  front  as  it  breaks  out  through  the  surface  of 
the  explosive  charge  in  the  modelling  results,  the  apparent  detonation  velocities  were 
determined  and  plotted  in  Figure  14,  together  with  the  experimental  results  [Dorsett 
and  Katselis  1999;  Dorsett  and  Jones  2001].  It  is  seen  that  the  hydrocode  model  can 
reproduce  accurately  the  experimental  results. 


Figure  14.  Comparison  of  predicted  and  measured  apparent  detonation  velocity  along  the  charge 
axis 

A  typical  pressure  contour  plot  generated  by  a  LS-DYNA  simulation  of  detonation  in 
water  is  shown  in  Figure  15,  also  showing  finite  element  meshes.  Figure  16  presents  the 
results  of  the  simulation  after  the  resolved  reaction  zone  detonation  in  the  PBXW- 
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115  (Aust)  has  rxin  approximately  32  ps.  The  calculated  contours  are  of  pressure  in  1.4 
GPa  intervals.  The  positions  of  the  measured  shock  front  and  the  bubble 
(explosive/ water)  interface  (after  correction  for  index-of-refraction  effects  [Craig  et  al 
1978])  are  included  as  open  circles  [Dorset!  and  Katselis  1999;  Dorset!  and  Jones  2001]. 
The  results  of  the  simulation  after  the  detonation  has  run  approximately  41  |as  and 
simulated  water  shock  profiles  off  the  end  of  the  charge  at  56  ps  together  with  the 
experimental  data  are  shown  in  Figures  17  and  18.  It  is  seen  that  the  simulated  shock 
and  explosive/ water  interface  are  in  excellent  agreement  with  the  experimental 
observation. 


The  good  agreement  between  the  hydrocode  modelling  and  experimental  observation 
showed  that  the  Ignition  and  Growth  Reactive  Model,  although  determined  from  the 
detonation  velocity  versus  diameter  effect  experimental  data  for  unconfined 
explosives,  could  faithfully  model  the  reaction  rate  characteristics  for  experiments 
carried  out  in  other  configurations.  This  gives  us  confidence  in  applying  the  mentioned 
models  to  simulation  of  the  mid-scale  underwater  tests  reported  by  Wilkinson  [2001] 
which  are  presented  in  the  following  section. 
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Time  a  31 39  |LS 


Figure  16.  The  simulation  results  by  LS-DYNA  (solid  lines)  of  the  PBXW-115(Aust) 
aquarium  test  (open  circles  represent  the  experimental  data) 


Figure  17.  Simulated  water  shock  profiles  after  the  resolved  reaction  zone  detonation  has  run 
approximately  41  jus  (open  circles  represent  the  experimental  data) 
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represent  the  experimental  data) 


3.6  Simulation  of  mid-scale  underwater  tests 

It  is  expected  that  the  Ignition  and  Growth  Reactive  Model  and  CPeX  Model  described 
above  using  data  from  small-scale  tests  can  be  applied  to  large-scale  devices.  To  test 
this  extrapolation,  both  the  Ignition  and  Growth  Model  in  LS-DYNA  and  the  CPeX 
Model  in  DYNA2D  [Kennedy  and  Jones  1993]  have  been  used  in  the  simulation  of  the 
mid-scale  underwater  tests  reported  by  Wilkinson  [2001].  Measurement  of  the 
underwater  explosive  performance  parameters  was  obtained  through  the  use  of  an 
array  of  pressure  gauges.  The  gauge  stations  were  placed  at  2m,  6m,  11m,  16.5m  and 
22m.  The  shock  wave  similitude  equation  parameters  for  peak  pressure  (Pm)  derived  by 
Wilkinson  [2001]  is  given  by 

p,„  =  50.42''-'^  (18) 

where  the  variable  Z  represents  the  scaled  distance  (R/W^/^),  R  (m)  is  the  distance  from 
the  centre  of  the  charge  and  W  (kg)  is  the  charge  mass.  Note  that  the  similitude 
parameters  are  explosive  specific  and  are  generally  specified  for  spherical  charges  with 
radius  of  R. 

Two  configurations  were  used  in  the  simulations.  The  first  was  used  to  study  the 
initiation  process  and  the  intermediate  expansion  of  the  water  shock  and  bubble  with 
both  the  Ignition  and  Growth  Model  and  the  CPeX  Model,  and  the  second  was  used  to 
examine  the  longer-term  effects  in  underwater  detonation  with  only  the  Ignition  and 
Growth  Model.  The  size  of  the  first  calculation  configuration  was  set  to  ^>530mm  x 
L1200mm  and  the  second  to  a>530mm  x  LSOOOmm  in  axisymmetric  cylindrical 
geometry.  The  dimensions  of  the  PE4  booster  and  PBXW-115  charge  are  34  and  126 
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mm  in  radius  and  68  and  400mm  in  length,  respectively.  The  thickness  of  the  fibreglass 
confinement  case  is  9mm.  Due  to  enormous  element  distortions  at  later  time,  the 
confinement  case  was  not  included  in  the  second  simulation  configuration.  The 
material  properties  and  equations  of  state  for  PE4,  PBXW-115(y4MsO  and  water  used  in 
the  calculation  are  the  same  as  those  used  in  the  aquarium  test  simulation.  The  model 
parameters  for  the  fibreglass  case  are  given  in  Table  8.  Figure  19  shows  the  mesh  of  the 
model  in  axisymmetric  geometry. 


Figure  19  Geometry  and  finite  element  of  the  model 


Table  8.  Material  and  Gruneisen  parameters  for  fibreglass 


Yield  stress 

0.2GPa 

Shear  modulus 

15GPa 

po  (g/cm3) 

1.70 

To 

1.01 

C  (mm/ps) 

3.016 

S 

1.005 

Figure  20  shows  the  initiation  process  and  the  intermediate  expansion  of  the  water 
shock  and  bubble.  It  should  be  noted  that  the  images  from  10  to  80  microseconds  are 
plotted  using  the  linear  pressure  scale  and  the  images  at  120  and  180  microseconds  are 
plotted  using  a  log  scale.  Good  agreement  is  seen  between  the  Ignition  and  Growth 
Model  and  the  CPeX  Model.  The  maximum  pressure  in  the  PBXW-115{>lusf)  occurs  at 
the  end  of  the  booster,  but  since  the  curvature  of  the  resulting  detonation  wave  is  too 
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large  for  it  to  be  sustained,  the  shock  attenuates  down  the  center  line  until  at  about  the 
center  of  the  charge,  the  shock  pressure  begins  to  build  up  again  as  its  curvature 
steadily  decreases.  The  detonation  propagates  most  strongly  in  a  45  degree  cone  out 
from  the  booster,  again  because  of  dependence  on  the  curvature. 


Figure  20.  Calculated  pressure  contours  at  different  times  with  I  &  G  model  plotted  on  the  left 
side  of  each  diagram  while  the  CPeX  model  is  plotted  on  the  right 

Modelling  the  far  field  late-time  effects  of  underwater  detonation  for  non-ideal 
explosives  such  as  PBXW-115  is  a  challenging  task.  To  this  end,  to  avoid  the  enormous 
element  distortion,  we  have  removed  the  inner  explosive  and  water  parts  when  the 
shock  wave  passed  this  portion.  Although  it  is  understood  that  doing  so  ignores  the 
momentum  behind  the  shock,  the  smooth  transition  curve  predicted  by  the  Ignition 
and  Growth  Model  in  LS-DYNA  as  shown  in  Figure  21  is  encouraging.  In  this  way,  the 
problem  can  rim  up  to  2.49m  from  the  charge  centre.  For  comparison,  the  peak 
pressure  fitted  to  equation  (18)  is  also  included  in  the  figure.  The  figure  shows  peak 
pressure  against  scaled  distance  Z=L/W^^,  where  L  (m),  is  the  distance  from  the  centre 
of  the  charge,  and  W  (kg)  is  the  charge  mass.  There  is  excellent  agreement  between  the 
LS-DYNA  simulation  and  the  experimental  data  from  the  aquarium  test  [Dorsett  and 
Jones  2001]  and  mid-scale  underwater  test  [Wilkinson  2001].  This  indicates  that  the 
Ignition  and  Growth  Model  can  capture  accurately  not  only  the  physics  of  detonation 
inside  the  explosives  but  also  the  flow  properties  in  the  region  between  the  outer 
boimdary  of  the  aquarium  and  a  scaled  distance  of  0.68.  It  is  worthy  of  note  that  there 
are  no  published  experimental  data  for  peak  pressures  in  this  region.  Figure  22  plots 
the  calculated  peak  pressures  versus  time  at  the  charge  edge  and  in  water.  The 
computed  peak  pressure  histories  at  various  locations  in  the  water  are  repeated  for 
clarity  in  Figure  23.  These  figures  show  the  gradual  decay  in  the  peak  pressure  with 
increasing  distance,  along  with  the  gradual  increase  in  arrival  and  duration  times. 
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Figure  21.  A  comparison  of  predicted  peak  overpressure  with  experimental  data 
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Figure  22  Peak  pressure  versus  time  at  the  charge  edge  and  in  water 
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Figure  23  Peak  pressure  histories  at  various  locations  within  water 

The  limitation  of  the  above  method  is  that  we  deleted  the  inner  portion  to  get  the 
problem  to  run  without  considering  the  perturbation  that  is  introduced  when  the 
elements  were  deleted.  In  that  way,  we  generated  a  rarefaction  that  will  degrade  the 
pressure  wave  from  behind,  limiting  how  far  it  can  run  before  being  totally  eroded. 

Currently  we  are  exploring  two  other  possible  approaches  to  tackle  the  problem.  It 
would  be  desirable  to  set  up  the  problem  from  the  start  with  three  parts,  one  for  the 
high  explosive  (HE),  one  for  the  inner  portion  of  the  water  that  will  be  deleted,  and  one 
for  the  outer  portion  of  water  that  will  remain. 

The  first  approach  is  to  run  the  simulation  first  and  record  the  pressure  time  history  in 
an  element  on  the  boundary  between  the  inner  and  outer  water.  If  the  HE  and  inner 
water  are  intended  to  be  deleted  at  time  ti,  the  first  simulation  needs  to  run  past  ti,  say 
to  t2  (t2  is  slightly  greater  than  ti).  An  exponential  function  then  needs  to  be  fitted  to  the 
pressure-time  data  between  ti  and  t2,  in  the  form 

P  =  Pq  (19) 

where  Po  and  A  are  the  fitting  parameters.  We  would  then  have  to  rerun  the  simulation 
setting  the  termination  time  to  ti,  though  now  with  a  ’^LOAD_SEGMENT_SET 
applying  a  pressure  load  curve  to  the  interface  segments  between  the  inner  and  outer 
water  elements,  with  the  birth  time  of  the  load  being  set  to  ti.  The  pressure-time  points 
for  the  load  curve  would  then  be  generated  using  the  function 

P  =  Po 
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When  this  second  simulation  stops  at  h,  we  would  then  restart  it  deleting  the  HE  and 
inner  water  parts.  The  pressure  boundary  will  automatically  turn  on  by  itself. 
Obviously,  the  artificial  pressure-time  load  is  not  going  to  be  exact  all  the  way  along 
the  new  boundary  (unless  the  shape  of  the  inner/ outer  water  interface  could  be 
matched  to  a  pressure  isobar  at  h),  but  the  disturbance  should  be  less  than  doing 
nothing  at  all.  If  there  is  a  particular  direction  where  it  is  more  important  to  get  the  far- 
field  correct,  the  pressure-time  history  should  be  recorded  there. 

A  second  way  to  approach  the  problem  is  to  switch  the  HE  and  inner  water  to  a  rigid 
body  at  time  ti  using  *DEFORMABLE_TO_RIGID_AUTOMATIC.  This  approach 
should  work  best  if  a  right  switch  time  can  be  picked  at  which  the  particle  velocities  in 
the  HE  and  the  inner  water  are  as  low  as  possible  (the  final  rigid  body  is  given  the 
centre-of-mass  velocity  of  the  parts  it  is  replacing). 

Work  on  investigating  the  above  two  possible  approaches  for  simulating  the  far  field 
for  late-time  effects  of  underwater  detonation  using  the  Ignition  and  Growth  Model  is 
currently  underway. 


4.  Conclusions  and  Future  Directions 

From  the  modelling  results  presented  using  kinetic  CHEETAH  and  the  simulation 
results  of  detonation  velocity  tests,  aquarium  tests  and  the  mid-scale  underwater  tests 
with  the  parameters  developed  for  both  the  Ignition  and  Growth  Reactive  Model  and 
CPeX  Model,  the  following  conclusions  have  been  reached. 

•  Kinetic  CHEETAH  can  predict  the  trend  of  the  detonation  velocity  versus 
diameter  effect.  The  square  root  pressure  dependence  provides  a  better  fit  to  the 
detonation  velocity  data. 

•  Preliminary  results  using  kinetic  CHEETAH  indicate  that  the  detonation 
velocity  and  the  critical  diameter  are  sensitive  to  the  assumed  AP 
decomposition  rate.  This  might  suggest  that  the  observed  differences  between 
the  US  and  the  Australian  composition  are  due  to  the  particle  size  of  the  AP 
rather  than  the  RDX.  Further  investigation  of  the  effects  of  the  assumed  AP 
decomposition  rate  on  the  detonation  velocity  as  a  function  of  charge  diameter 
will  be  carried  out  in  the  future. 

•  The  Igmtion  and  Growth  Reactive  Model  can  capture  the  physics  of  detonation 
inside  the  explosives  accurately,  and  successfully  simulate  successive  positions 
of  bubble  and  shock  front  resulting  from  a  detonation  wave  propagating  down 
the  explosive  cylinder  observed  in  the  aquarium  tests. 

•  Both  the  Ignition  and  Growth  Reactive  Model  and  the  CPeX  model,  although 
determined  from  the  detonation  velocity  versus  diameter  effect  experimental 
data  for  unconfined  explosives,  could  faithfully  model  the  reaction  rate 
characteristics  for  experiments  carried  out  in  other  configurations.  This  has 
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been  validated  by  comparing  the  predictions  using  both  models  against 
experimental  data  for  detonation  of  charges  confined  in  2.5  mm  and  3  mm  thick 
brass. 

•  The  good  agreement  between  the  Ignition  and  Growth  Model  and  the  CPeX 
Model  for  the  simulation  of  mid-scale  underwater  tests  implies  that  these  two 
models  using  parameters  derived  from  small-scale  tests  could  be  applied  to 
large-scale  devices.  We  have  used  a  tentative  method  to  model  the  far  field  late¬ 
time  effects  of  underwater  detonation  by  deleting  the  inner  explosive  and  water 
parts  when  the  shock  wave  passed  the  portion.  Although  it  is  understood  that 
the  method  ignores  the  momentum  behind  the  shock  without  considering  the 
perturbation  that  is  introduced  when  the  elements  were  deleted,  the  smooth 
transition  curve  predicted  between  the  model  and  the  available  experimental 
peak  overpressure  data  is  promising.  We  have  also  outlined  two  other  possible 
approaches.  Work  on  investigating  these  two  possible  approaches  for 
simulating  the  farther  field  for  late-time  effects  of  imderwater  detonation  using 
the  Ignition  and  Growth  Model  is  currently  underway. 

•  The  LS-DYNA  simulations  of  shock  front  curvature  for  both  PBXN-111  and 
PBXW-115{Ausf)  unconfined  charges  using  Ignition  and  Growth  Reactive 
Model  are  presented,  where  the  simulated  results  for  PBXN-111  are  compared 
with  the  available  experimental  data.  The  average  sonic  reaction  zone  lengths 
for  both  PBXN-111  and  PBXW-115(/lMsf)  unconfmed  charges  are  calculated 
based  on  the  simple  wave  curvature  theory  by  Souers  and  Garza.  Given  that 
there  is  the  lack  of  experimental  wave  curvature  data  for  PBXW-115(Aust),  the 
simulation  results  presented  herein  could  provide  some  insight  into  the 
understanding  of  the  detonation  process  for  this  non-ideal  explosive. 

•  As  part  of  the  program  of  evaluating  the  ultrafine  aluminium  powder  'Alex'  in 
explosive  formulations,  aquarium  tests  of  Alex-based  PBXW-115(Awsf)  will  be 
carried  out  to  measure  Alex  reaction  rates  in  polymer-bonded  explosives.  These 
experimental  results,  together  with  the  results  of  VoD  and  plate  dent  depth 
tests  and  aquarium  test  of  Tritonal  (80:20  TNT/Al)  explosive  formulations 
containing  'Alex'  [Cliff  et  al  2002],  will  be  used  to  calibrate  the  Ignition  and 
Growth  model  in  LS-DYNA  to  study  the  role  of  aluminium  and  particle  size 
effects. 
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perchlorate,  25  %  aluminium,  20%  RDX  and  12%  HTPB  binder,  has  been  studied  using  Kinetic  CHEETAH  (the 
Lawrence  Livermore  National  Laboratory  CHEETAH  2.0  code),  the  three-term  "Ignition  and  Growth"  Model 
incorporated  into  the  explicit  finite  element  hydrocode  LSTC-DYNA,  and  the  CPeX  Reactive  Model  based  on 
the  small  divergent  detonation  theory  of  Wood  and  Kirkwood  in  DYNA2D.  This  report  firstly  focuses  on  a 
series  of  simulations  performed  using  Kinetic  CHEETAH  based  on  the  Wood-Kirkwood  detonation  theory 
using  a  pressure  exponent  of  2  in  the  pressure-dependent  rate  law.  It  was  found  that  CHEETAH  could  predict 
the  trend  of  the  detonation  velocity  as  a  function  of  charge  diameter,  but  overestimated  the  detonation 
velocities.  The  agreement  was  improved  with  further  parameter  adjustment,  by  decreasing  the  pressure 
exponent  in  the  rate  law  from  2.0  to  1.0,  and  then  0.5.  The  reaction  zone  widths  over  a  wide  range  of  charge 
radius  were  also  computed.  Attention  was  then  turned  to  hydrocode  modelling,  with  particular  interest  in 
developing  reactive  models  for  PBXN-111  and  PBXW-115  (Aust).  Both  the  "Ignition  and  Growth"  Model  and 
the  CPeX  Model  were  calibrated  against  the  experimentally  observed  dependence  of  detonation  velocity  on 
charge  diameter  in  unconfined  charges  of  both  PBXW-115  (Aust)  and  US  PBXN-111.  These  two  reactive  models 
were  validated  by  comparing  their  predictions  against  experimental  data  for  detonation  of  charges  confined  in 
2.5  mm  and  3  mm  thick  brass.  The  Ignition  and  Growth  Model  was  then  applied  successfully  to  the  simulation 
of  aquarium  tests.  Finally,  to  test  whether  the  two  models  using  parameters  derived  from  small-scale  tests  can 
be  applied  to  large-scale  devices,  these  two  models  were  applied  to  the  simulation  of  mid-scale  imderwater 
tests  of  PBXW-115  {Aust),  and  comparisons  with  available  data  are  made.  The  LS-DYNA  simulations  of  shock 
front  curvature  for  unconfined  charges  are  also  presented  and  the  relationships  between  the  radii  of  curvature 
for  the  detonation  fronts  and  self-propagating  detonation  velocities  are  discussed. 
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